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Meet  the 

Cluster  Orbits  With  Perturbations  Of  Keplerian  Elements  (COWPOKE) 

Equations 

Chris  Sabol*,  Craig  A.  McLaughlin^  and  K.  Kim  Luu* 

Recent  developments  have  indicated  that  it  is  possible  to  express  the  relative  equations  of 
motion  for  space  objects  in  non-circular  orbits  using  mean  Keplerian  elements  and  low  order 
expansions.  This  paper  provides  the  initial  derivation  of  one  such  effort  known  as  the  Cluster 
Orbits  With  Perturbations  Of  Keplerian  Elements  (COWPOKE)  equations.  Given  mean 
Keplerian  elements  and  element  differences,  the  COWPOKE  equations  describe  spherical 
radial,  cross-track,  and  along-track  separations  of  the  satellites  as  an  explicit  function  of  time. 

The  framework  of  the  equations  allows  for  very  high  eccentricity  reference  orbits  and  for  the 
inclusion  of  dynamic  perturbations.  Test  cases  using  two-body  dynamics  show  the  utility  of 
tiiis  approach. 

INTRODUCTION 

Cluster  orbits  are  defined  as  the  relative  trajectories  of  objects  traveling  through  space  in  close  proximity  to  each 
other.  Clusters  occur  when  satellites  in  similar  orbits  approach  each  other,  groups  of  satellites  fly  in  formation  to 
perform  specific  mission  functions,  or  several  objects  are  launched  or  deployed  from  the  same  object.  Each  case  has 
its  own  unique  nuances,  but  the  dynamics  of  how  the  objects  move  relative  to  each  other  are  functionally  the  same. 

The  number  of  satellites  in  the  geosynchronous  belt  has  been  steadily  increasing  for  almost  40  years.  While  the 
number  of  satellites  has  been  increasing,  the  space  in  the  geosynchronous  belt  has  not.  This  has  led  to  more  clusters 
of  satellites  operating  in  close  proximity  of  each  other  and  often  unintentionally  passing  within  kilometers  of  ftieir 
neighbors*.  Clusters  of  objects  at  this  great  altitude  have  been  a  long-standing  challenge  to  space  surveillance  where 
the  satellites  can  be  mistaken  for  each  other. 

In  recent  years,  there  has  been  increasing  interest  in  the  use  of  satellites  flying  in  formation.  Several  missions  and 
mission  statements  have  identified  formation  flying  as  a  means  of  reducing  cost  and  adding  flexibility  to  space  based 
programs  or  to  accomplish  goals  that  are  not  possible  or  very  difficult  to  accon:q)lish  with  a  single  satellite.  These 
missions  include  NASA's  Earth  Observing-1  flying  in  formation  wifti  Landsat-7  and  several  European  missions.  In 
addition  NASA  has  developed  over  20  concepts  for  future  missions  involving  formation  flying.  Many  of  these 
missions  involve  highly  eccentric  orbits. 

Dining  ballistic  missile  launches,  several  objects  can  achieve  low  earth  orbit  including  the  reentry  vehicle,  booster, 
fairings^  and  the  possibility  of  several  balloon-like  decoys  attenqjting  to  confuse  missile  defense  systems.  It  is  the 
responsibility  of  the  missile  defense  system  to  track  the  cluster  of  objects  and  quickly  identify  the  reentry  vehicle. 

The  Cluster  Orbits  with  Perturbations  of  Keplerian  Elements  (COWPOKE)  equations  can  provide  the  theoretical 
foundation  for  analysis  tools  supporting  a  variety  of  applications.  Optical  space  surveillance  often  reveals  multiple 
satellites  in  a  single  telescope  field  of  view.  An  understanding  of  the  relative  dynamics  may  allow  for  satellite 
identification  based  solely  on  relative  position.  Further  research  may  apply  these  results  to  relative  orbit 
determination  of  objects  in  a  cluster  or  during  close  approach  encounters.  This  approach  may  prove  to  be  more 
accurate  than  performing  orbit  determination  on  the  individual  objects  and  determining  differences.  The  equations 
might  also  be  used  to  support  the  missile  tracking  experiments.  If  tracking  sensors  have  limited  fields  of  view,  they 
may  be  required  to  track  ballistic  clusters  individually;  meaningful  relative  equations  of  motion  may  simplify  the 
transition  from  object  to  object.  Even  if  objects  are  tracked  simultaneously,  the  relative  equations  of  motion  can  be 
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used  as  part  of  the  discrimination  process  to  identify  balloon-like  decoys.  The  COWPOKE  equations  would  also  be 
extremely  relevant  to  satellite  formation  flying  worl^  such  as  the  initial  Air  Force  Research  Laboratory’s  TechSat  21 
program,  for  formation  flying  design,  analysis,  and  guidance  and  control  applications.  While  a  great  deal  of  work 
has  already  been  completed  in  support  of  formation  flying  missions,  none  provide  the  level  of  intuition  and 
understanding  as  one  that  uses  Keplerian  elements.  In  addition,  none  of  the  work  is  accurate  for  long-term  relative 
motion  for  orbits  with  eccentricities  up  to  0.7. 

A  simple  set  of  equations  describing  the  relative  motion  of  spacecraft  clusters  is  needed  for  analysis,  design,  and/or 
tracking  of  clusters.  Much  of  die  previous  work  in  this  field  has  relied  on  using  Hill's  equations^  (also  known  as  the 
Clohessy-Wiltshire  equations^).  Hill’s  equations  describe  the  relative  motion  of  spacecraft  using  a  spacecraft- 
centered  coordinate  system.  However,  Hill's  equations  assume  that  the  reference  orbit  is  circular,  the  objects  are 
close  together,  and  there  are  no  perturbations  to  simple  two-body  motion.  Several  researchers  have  pointed  out  the 
severe  limitations  inherit  in  these  assunqjtions.  To  rectify  this,  Gim  and  Alfnend^  used  energy  methods  to  develop  a 
state  transition  matrix  to  describe  relative  motion  in  non-circular  orbits  under  the  influence  of  perturbations.  While 
closed  form  analytic  solutions  in  a  transformed  variable  space  are  valuable  for  analysis,  these  approaches  may  be 
awkward  for  many  applications  since  transformations  are  required  to  go  from  the  canonical  variable  space  to  more 
traditional  representations  of  satellite  orbits  such  as  Keplerian  elements.  Additionally,  most  working  level  engineers 
lack  the  backgroxmd  to  inclement  this  approach.  Previous  work  by  Garrison  et  al.^  developed  equations  of  motion 
for  elliptical  orbits  in  terms  of  the  true  anomaly  instead  of  time.  Also,  these  equations  were  developed  with 
rendezvous  in  mind  and  are  developed  only  to  second  order  in  eccentricity.  In  addition,  Melton  developed  a  state 
transition  matrix  for  relative  motion  in  eccentric  orbits  that  is  time  dependent.  Melton  relied  on  some  of  the  same 
foimdations  that  are  used  in  this  paper,  but  the  development  was  only  to  second  order  in  eccentricity  and  is  not 
accurate  for  orbits  of  high  eccentricity.  Baoyin  et  al,*^  did  similar  work  but  assumed  matching  orbital  periods;  this 
work  also  provides  physical  insight  into  formation  flying  in  near  circular  orbits.  This  paper  develops  physically 
meaningful  equations  of  relative  motion  for  space  objects  in  non-circular  orbits  using  Keplerian  elements. 

Previous  work  by  the  authors’  led  to  the  development  of  a  set  of  equations  that  describe  the  first  order  effects  of 
Earth  oblateness  on  the  relative  motion  of  objects  in  circular  orbits.  In  addition,  fliey  developed  a  single  set  of 
equations  to  describe  the  effects  of  Earth  oblateness  for  polar  orbits.  Then,  they  further  developed  the  equations  to 
describe  the  motion  for  all  inclinations*.  Finally,  they  examined  the  long-term  evolution  of  the  relative  motion  for 
cireular  orbits  and  presented  an  approach  to  describe  relative  motion  of  satellites  in  elliptieal  orbits  without 
perturbations  and  assuming  matching  periods®.  This  last  step  provided  the  building  blocks  for  the  COWPOKE 
equations. 

In  the  past,  equations  of  relative  motion  such  as  Hill’s  equations  were  developed  by  differencing  the  quasi-inertial 
equations  of  motion  and  mapping  those  differences  into  the  rotating  coordinate  system.  The  result  was  a  set  of 
equations  that  describe  the  relative  differences  between  the  satellites  typically  in  terms  of  radial,  cross-track,  and 
along-track  coirqjonents.  Rather  than  using  the  traditional  algebraic  approach,  the  COWPOKE  equations  are 
developed  using  a  geometric  approach.  Here,  the  geometric  properties  and  definitions  of  Keplerian  orbital  elements 
will  be  used  directly  to  map  orbital  element  differences  into  the  radial,  cross-track,  and  along-track  relative  motion. 
In  addition  to  the  increased  intuitiveness  of  a  geometric  approach,  by  using  Keplerian  orbital  elements  and 
differences  in  those  orbital  elements,  we  can  take  advantage  of  existing  perturbation  model  development  and 
incorporate  meaningful  dynamics  into  the  relative  equations  of  motion  much  easier  than  with  the  algebraic 
approaches. 

There  are  challenges  to  this  approach.  First  and  foremost,  die  geometric  properties  of  Keplerian  elements  are  a 
function  of  the  tme  anomaly,  wWch  has  non-uniform  variation  with  time.  Ideally,  analytical  solutions  are  expressed 
in  terms  of  the  mean  anomaly,  which  does  vary  linearly  with  time  (outside  the  influence  of  perturbations).  The 
relationship  between  the  mean  anomaly  and  eccentric  anomaly  is  given  by  the  well-  known  transcendental  Kepler’s 
equation.  This  research  uses  a  series  expansion  to  approximate  the  trae  anomaly  as  a  function  of  time.  A  relation  of 
this  sort  was  successfully  employed  as  part  of  the  building  blocks  of  the  COWPOKE  equations. 

Using  that  one  relationship,  one  can  constract  the  relative  motion  of  formation  flying  satellites  by  reconstmcting  the 
orbits  individually  based  on  their  Keplerian  elements  and  then  differencing  the  two.  This  paper  develops  true 
equations  of  relative  motion  for  clusters  of  objects.  And,  unlike  the  previous  work,  this  research  incorporates 


perturbation  models  and  be  generalized  to  account  for  small  differences  in  semimajor  axis.  This  results  in  a  set  of 
equations  that  describe  the  radial,  cross-track,  and  along-track  differences  between  two  space  objects  in  formation  or 
a  cluster  based  on  Keplerian  orbital  element  differences. 

APPROACH 

The  formulation  of  the  COWPOKE  equations  has  three  major  challenges.  The  first  is  choosing  a  suitable  reference 
fi-ame  for  the  relative  equations  of  motion  and  expressing  that  frame  in  terms  of  Keplerian  elements  and  element 
differences.  The  second  is  representing  the  geometrical  hue  anomaly  as  a  function  of  time.  The  final  challenge  is 
incorporating  relevant  perturbation  effects.  Only  the  first  two  challenges  are  addressed  in  this  work. 

Classical  relative  motion  approaches  use  a  Cartesian  reference  firame  to  describe  relative  motion.  The  conqjonents  of 
this  fi'ame  are  radial,  cross-track,  and  along-track  differences.  A  limitation  of  this  approach  is  that  significant  cross¬ 
mapping  between  components  occurs  when  the  separation  distance  between  the  satellites  is  not  small.  For  instance, 
if  two  satellites  travel  together  in  a  circular  orbit  but  are  separated  by  0.01  rad  in  true  anomaly,  the  relative  position 
of  one  satellite  with  respect  to  the  other  will  have  an  along-track  separation  as  anticipated  but  will  also  have  a  radial 
conqjonent  despite  the  fact  that  die  satellites  are  at  the  same  altitude  and  within  the  same  orbit.  Relative  position 
difference  components  that  are  artifacts  of  die  chosen  coordinate  fi:ame,  such  as  the  one  outlined  in  the  preceding 
exan^le,  can  greatly  conqilicate  the  analysis  of  the  relative  motion. 

A  second  coordinate  firame  choice  is  to  use  spherical  separations.  Here,  the  relative  position  difference  is  described 
by  an  altitude  difference  firom  a  sphere  having  the  radius  of  the  reference  satellite  and  angular  components 
perpendicular  to  and  along  the  reference  satellite’s  direction  of  motion  projected  onto  the  sphere.  These 
conqionents,  dr.  Set,  and  Sat,  are  illustrated  in  Figure  1.  The  spherical  reference  fi-ame  is  well-suited  to  describing 
the  position  of  satellites  in  circular  or  near-circular  orbits;  however,  when  die  orbital  eccentricity  becomes  large,  the 
spherical  reference  frame  suffers  from  the  same  limitations  as  the  Cartesian  system. 

A  third  option  to  describe  the  relative  motion  of  satellites  in  elliptic  orbits  is  to  use  an  ellipsoidal  reference  system. 
This  is  similar  to  the  spherical  system  except  that  the  components  are  mapped  along  an  ellipsoid  rather  than  a  sphere. 
This  method  has  geometric  challenges  describing  altitude  variations  and  can  still  result  in  difficult  to  understand 
position  differences. 

For  this  formulation,  we  chose  the  spherical  reference  frame  since  it’s  slightly  better  than  Cartesian  for  describing  the 
orbital  motion  and  provides  the  physical  insight  we  desire.  Additionally,  it  is  a  simple  geometric  mapping  between 
spherical  and  Cartesian  given  the  quantities  described  in  this  work  if  Cartesian  coordinates  are  desired.  In  fact,  the 
desired  reference  firame  will  likely  be  a  function  of  application,  but  it  is  hoped  that  enough  information  is  provided 
here  so  these  equations  can  be  easily  reformulated  in  other  reference  frames. 


Figure  1:  The  Spherical  Reference  Frame  for  Relative  Motion 


The  spherical  coordinate  system  used  here  can  be  realized  using  position  and  velocity  vectors  of  the  two  satellites: 
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Here,  the  auxiliary  vector  quantities,  denoted  by  primes,  are  introduced  so  right  planar  triangles  can  be  used  in  the 
determination  of  desired  angular  values.  The  spherical  components  can  be  expressed  in  terms  of  Keplerian  elements 
by  using  the  following  definitions**’: 
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However,  the  resulting  equations  become  difficult  to  manipulate.  Instead,  a  simple  geometrical  mapping  of 
Keplerian  element  and  element  differences  into  the  spherical  components  of  relative  motion  can  be  achieved  with  the 
following: 

_  (a,+Sa)(l-(e,+Sef)  a,(l-e,^) 

l  +  (e^  +  Se)cos(f^  +  Sf)  l  +  eiCOs(/) 

=P  =  -^sin(t,)cos(0,  +  +  /  +  ^)  +  sin(6)j  +  +  /  +  c^)  (3) 


=a  =  {So}  + df)cos{8i)  +  dD.cos{i^) 

/  ^ 

The  spherical  cross-track  and  along-track  terms  are,  in  fact,  approximations  accurate  to  first  order  in  the  Keplerian 
elements  differences.  Figure  2  illustrates  how  these  terms  are  derived. 


Figure  2:  Spherical  Conqionents  in  Terms  of  Keplerian  Elements  and  Element  Differences 


We  can  further  simplify  the  radial  expression  to  first  order  in  the  Keplerian  element  differences: 

l-fe,cos(/i)  (lH-ejCos(/,))^  (1-i-e,  cos(/,))^ 

From  these  Keplerian  element  representations  of  the  relative  motion,  we  can  derive  some  physical  insight  into  how 
orbital  differences  result  in  relative  motion.  If  we  assume  a  circular  reference  orbit,  Eqs.  (3)  and  (4)  resemble  the 
solutions  to  Hill’s  equations  where  the  cross-track  motion  becomes  a  single  oscillation,  the  altitude  difference  maps 
directly  into  the  radial  motion  with  a  periodic  variation  introduced  by  the  eccentricity,  and  the  altitude  difference 
couples  into  the  along-track  conqjonent  through  the  true  anomaly  resulting  in  a  combination  of  secular  and  periodic 
effects’. 

While  the  equations  above  are  useful  and  do  help  provide  physical  insight  into  the  relative  motion,  they  are 
expressed  in  terms  of  true  anomaly.  To  express  the  relative  motion  as  a  function  of  time,  one  must  express  the  above 
equations  in  terms  of  mean  anomaly.  At  a  fundamental  level,  this  entails  finding  a  suitable  approximate  solution  to 
Kepler’s  equation.  Battin  provides  several  methods  for  doing  so*'*.  For  the  initial  COWPOKE  development,  it  was 
decided  to  use  a  Fourier-Bessel  expansion  of  die  true  anomaly  in  terms  of  the  mean  anomaly  and  eccentricity. 

Using  the  geometric  properties  of  the  true  and  eccentric  anomalies  and  Kepler’s  equation,  the  following  relationship 
can  be  found: 
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This  can  be  expressed  in  terms  of  a  Fourier  cosine  series  that  introduces  the  Bessel  functions, 


=  Vl-e'  +  £5*  cos(W) 

dM 


\  *=i 

2 


VW 


.  -80  = 


Vl-C^ 


(6) 


j  ( /;c~)"y  (~^r(~^) _ 


2w 


i-vr 


.A'*"' = 


/  \|*+n| 

'  e  1 


l^  +  n 
/! 


]^(jA:  +  w|+2/-y) 


M 


Integrating  this  equation  with  respect  to  the  mean  anomaly  provides  an  expression  for  the  true  anomaly  in  terms  of 
the  mean  anomaly: 


00  1 

f=M+2Z- 

k=\ 


Ej.ete)/*'' 


sm{kM) 
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The  result  of  this  expansion  is  a  sine  series  with  coefficients  that  are  power  series  in  eccentricity.  The  lowest  order 
of  eccentricity  for  each  coefficient  series  is  k;  thus  the  upper  limit  for  k  can  be  chosen  based  on  the  eccentricity  of 
the  reference  orbit  and  the  desired  accuracy.  Figure  3  plots  eccentricity  to  the  (k+1)  power.  The  figure  is  meant  to 
provide  an  indicator  for  the  relative 'accuracy  of  Eq.  (7)  for  a  given  eccentricity  and  choice  of  k.  For  instance,  if  one 
had  a  references  eccentricity  of  0.3  and  desired  series  truncation  errors  below  1%,  then  one  would  choose  to  truncate 
the  series  at  k=4  or  the  fourth  power  of  eccentricity.  When  k=8,  one  sees  errors  below  10%  up  to  e=0.7 . 
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Figine  3:  Relative  Error  Order  of  Magnitude  for  Eccentricity  Series 


Evaluating  Eq.  (7)  up  to  the  ninth  frequency  term  (8M)  yields: 
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In  examining  Eq.  (8),  one  can  see  that  the  leading  term  in  each  of  the  eccentricity  series  does  not  decrease  for  all 
eccentricities;  however,  the  second  term  in  the  eccentricity  series  is  always  negative  and  serves  to  reduce  the 
magnitude  of  that  frequency’s  contribution  to  the  expansion  as  a  whole.  Thus,  the  anqrlitude  of  each  frequency  term 
decreases  as  the  frequency  increases  allowing  convergence  of  the  entire  expansion  if  enough  terms  are  included  in 
each  individual  eccentricity  series.  This  means  that  for  higher  values  of  eccentricity,  it  is  very  in^ortant  to  include  at 
least  two  or  more  terms  for  each  eccentricity  series.  In  Eq.  (8),  if  the  e’  and  e‘“  terms  were  left  out  of  the  IM  and  8M 
frequency  terms,  the  error  due  to  this  truncation  could  be  larger  than  if  the  IM  and  8M  terms  were  not  included  at  all. 
For  higher  values  of  eccentricity,  care  must  be  taken  to  make  certain  that  enough  terms  are  included  in  each 
eccentricity  series  to  ensure  one  is  not  adding  error  with  a  given  frequency  term.  The  good  news  is  that  each 
eccentricity  series  converges  fairly  quickly  so  that  powers  of  eccentricity  at  lower  frequencies  are  much  less 
important  than  at  higher  frequencies;  this  means  that  one  can  get  away  with  fewer  terms  in  each  eccentricity  series 
even  if  higher  frequency  terms  are  required.  Additionally,  for  smaller  values  of  eccentricity,  0.3  and  below, 
truncation  becomes  much  less  of  an  issue  since  the  e^  factor  greatly  reduces  the  impact  of  the  higher  order  terms  in 
each  eccentricity  series  and  each  higher  frequency. 


Truncation  issues  aside,  the  Fourier-Bessel  series  expansion  of  the  true  anomaly  is  an  important  component  to  the 
formulation  of  the  COWPOKE  equations;  however,  an  additional  step  must  be  taken  to  map  mean  anomaly 
differences  into  true  anomaly  differences.  This  is  acconqjlished  by  simply  applying  a  first  order  perturbation  to  Eq. 
(8): 
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The  general  formulation  includes  higher  order  element  difference  terms  and  coupling  terms  between  the  element 
differences  but  for  cluster  orbits,  including  the  first  order  terms  should  provide  accuracy  to  two  orders  of  magnitude 
smaller  than  the  separation  distance.  Including  only  the  first  order  difference  terms  for  the  true  anomaly  expression 
yields: 
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The  truncation  issues  for  Eq.  (10)  are  similar  to  those  of  Eq.  (8)  but  are  furflier  confounded  by  the  exponential 
factors  produced  when  taking  the  partial  derivative  with  respect  to  eccentricity.  Note  that  the  e’  coefficient  is  twice 
as  large  as  the  e’  coefficient  in  the  8Af  frequency  series  for  the  8e  conf  onent. 

For  unperturbed  satellite  motion,  we  can  express  the  mean  anomaly  as  a  function  of  time  SM  as  a  function  of  SM  at 
epoch,  semimajor  axis  difference,  and  time: 


If  one  chooses,  the  mean  motion  tenns  in  the  above  expression  can  be  replaced  by  a  second  order  perturbation  of  the 
mean  motion  with  respect  to  Sa  without  significant  error  for  most  applications: 
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Eqs.  (8),  (10),  and  (11)  can  now  be  substituted  into  Eq.  (3)  to  complete  the  COWPOKE  equations  for  unperturbed 
satellite  motion.  Here  are  the  COWPOKE  equations  with  only  first  order  eccentricity  terms  included  in  the 
expansion  of  the  true  anomaly  and  true  anomaly  difference: 
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Note  that  we  have  taken  the  linearized  form  of  the  radial  con:5)onent  as  presented  in  Eq.  (4).  From  Figure  3,  we 
would  expect  that  these  equations  are  accurate  to  the  1%  level  up  to  eccentricities  of  0.1. 


While  this  work  does  not  currently  include  perturbations  to  the  satellite  dynamics,  one  of  the  key  focuses  of  the 
COWPOKE  development  was  to  allow  for  the  inclusion  of  additional  force  model  effects.  This  can  be  acconqilished 
by  expressing  the  Keplerian  elements  and  element  differences  as  fimctions  of  time.  This  will  be  the  focus  of  future 
work. 


RESULTS 

Simulations  were  performed  to  quantify  the  error  inherent  in  the  approximations  used  in  the  COWPOKE  formulation 
for  two  test  cases.  The  first  test  case  employed  a  near-circular,  low-Earth  orbit  (e=0.01)  and  the  second  a  high 
altitude  eccentric  orbit  (e=0.7).  The  test  cases  were  performed  in  the  MATLAB  environment  using  two-body 
dynamics;  it  is  understood  that  that  these  simulations  only  measure  the  effectiveness  of  the  COWPOKE  equations  to 
model  two-body  motion  and  are  not  indicative  of  the  performance  for  real  satellite  motion. 

Truth  data  for  the  simulations  were  determined  by  calculating  the  mean  anomaly  of  each  satellite  as  a  function  of 
time,  converting  the  Keplerian  elements  to  position  and  velocity  vectors  (with  an  algorithm  taken  fi:om  Vallado"), 


and  then  mapping  the  position  differences  into  the  spherical  radial,  cross-track,  and  along-track  components 
described  in  Eq.  (1).  The  COWPOKE  results  come  from  a  direct  mapping  of  the  Keplerian  element  and  element 
differences  into  those  conqjonents  using  the  methods  described  in  the  previous  section.  It  should  be  noted  that  the 
angular  cross-track  and  along-track  values,  a  and  P,  were  used  for  conqiarison  purposes  in  the  simulations  rather 
than  the  arc-lengths,  dxt  and  Sat.  Some  additional  error  in  the  arc-lengths  will  be  present  due  to  die  error  in 
estimating  the  magnitude  of  the  position  vector  of  the  reference  satellite,  but  that  is  typically  the  same  order  of 
magnitude  or  smaller  than  the  error  in  the  angular  separation. 

Table  1  contains  the  orbital  elements  and  element  differences  for  the  near-circular,  low-Earth  orbit  (LEO)  case;  the 
reference  satellite,  satellite  1,  has  the  initial  Keplerian  elements  given  in  the  “Reference  Elements”  column  while  the 
second  satellite,  satellite  2,  has  the  initial  Keplerian  elements  of  the  “Reference  Elements”  plus  the  “Element 
Differences.”  The  simulation  span  covers  2  hours  or  just  over  one  orbital  period.  The  COWPOKE  formulation 
includes  only  first  order  eccentricity  terms  and  uses  fire  radial  conponent  approximation;  these  equations  were 
identical  to  Eq.  (13). 


Table  1:  Keplerian  Elements  and  Element  Differences  for  LEO  Test  Case 


Reference 

Element 

Elements 

Differences 

a 

7000  km 

0.01  km 

e 

0.01 

0.01 

i 

0.785  rad  (45  deg) 

O.Olrad 

Q 

0  rad  (0  deg) 

0.01  rad 

0) 

4.712  rad  (270  deg) 

0.01  rad 

Mo 

1.751  rad  (90  deg) 

0.01  rad 

Figures  4-6  plot  the  spherical  radial  difference,  angular  cross-track,  and  angular  along-track  separations, 
respectively.  Each  figure  actually  contains  two  plots:  the  first  shows  the  satellite  separations  described  by  the  trath 
orbits  and  the  COWPOKE  approximation,  and  the  second  shows  the  differences  between  the  truth  and  COWPOKE 
methods.  Thus,  the  second  plot  in  each  figure  is  the  error  inherent  in  the  given  COWPOKE  formulation. 

Figiue  4  shows  the  COWPOKE  radial  error  to  be  around  the  2%  level.  Figures  5  and  6  show  the  COWPOKE  cross¬ 
track  and  along-track  errors  are  below  1%.  It  is  interesting  to  note  that  the  radial  and  cross-track  conponent  errors 
have  frequencies  that  appear  to  be  once  per  orbit  while  the  along-track  component  has  an  error  frequency  of  twice 
per  orbit.  We  test  the  COWPOKE  approximations  directly  to  determine  the  dominant  sources  of  these  errors. 

Before  continuing  with  the  error  analysis,  it  is  useful  to  examine  Figures  4-6  and  observe  the  physical  insight  into  the 
relative  motion  we  can  gain  from  the  COWPOKE  equations.  Eq.  (4)  tells  us  fliat  the  relative  motion  in  the  radial 
direction  is  dominated  by  a  once  per  orbit  signature  due  to  the  eccentricity  differences;  the  contributions  from  the 
semimajor  axis  difference  are  small,  and  the  contribution  from  the  true  anomaly  difference  is  an  order  of  eccentricity 
smaller  than  the  contribution  from  the  eccentricity  difference.  Eq.  (3)  shows  that  the  cross-track  motion  is  conprised 
of  once  per  orbit  signatures  from  the  inclination  and  right  ascension  of  the  ascending  node  differences.  Eq.  (13) 
indicates  that  the  along-track  motion  is  offset  from  zero  due  to  the  right  ascension,  argument  of  perigee,  and  mean 
anomaly  differences  but  also  has  a  once  per  orbit  periodic  signature  with  an  anplitude  of  twice  the  reference 
eccentricity. 


Figure  6:  LEO  Test  Case  Angular  Along-Track  Differences  and  COWPOKE  Errors 


Figure  7  plots  the  error  in  the  first  order  Fourier-Bessel  series  expansion  of  the  true  anomaly  of  the  reference 
satellite.  From  Figure  3,  we  expect  to  see  around  1%  error,  and  we  do.  Also  note  the  twice  per  orbit  signature  in  the 
error  plot;  this  is  also  expected  since  the  expansion  only  include  once  per  orbit  terms.  Figure  8  plots  the  error  in  the 
approximation  of  the  true  anomaly  difference;  the  error  in  this  approximation  is  similar  to  the  error  in  the  true 
anomaly  approximation.  One  sees  that  the  dominant  part  of  the  along-track  error  is  due  to  the  COWPOKE 
approximation  error  of  the  true  anomaly  difference.  If  one  chose,  this  error  could  be  reduced  by  adding  additional 
terms  to  the  series  approximation  for  this  variable. 
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Figure  7:  True  Anomaly  Approximation  Error  for  LEO  Test  Case 
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Figure  8:  True  Anomaly  Difference  Approximation  Error  for  LEO  Test  Case 


Figure  9  plots  the  errors  in  the  COWPOKE  coordinate  system  approximations  for  the  LEO  test  case.  One  can  see 
that  the  dominant  part  of  the  radial  and  cross-track  errors  is  due  to  the  coordinate  system  approximations  made  in  the 
COWPOKE  formulations.  The  radial  conqjonent  error  can  be  reduced  by  using  the  analytical  formulation  for  the 
radial  difference  while  an  improved  description  of  the  cross-track  component  would  have  to  be  derived  to  reduce  the 
cross-track  error. 


Table  2  contains  the  orbital  elements  and  element  differences  for  the  highly  eccentric  Earth  orbit  (HEO)  test  case. 
The  simulation  spans  12  hours  or  just  over  one  orbital  period.  The  COWPOKE  formulation  includes  eighth  order 
eccentricity  up  to  frequency  6M  and  tenth  order  eccentricity  in  frequency  IM  and  8M;  these  are  all  of  the  terms 
included  in  Eqs.  (8)  and  (10).  This  test  case  does  not  use  the  radial  component  approximation. 


Reference 

Element 

Elements 

Differences 

a 

27000  km 

0.01  km 

e 

0.7 

0.01 

i 

0.785  rad  (45  deg) 

O.Olrad 

Q 

0  rad  (0  deg) 

0.01  rad 

(0 

4.712  rad  (270  deg) 

0.01  rad 

Mo 

1.751  rad  (90  deg) 

0.01  rad 

Figures  10-12  plot  the  spherical  radial  difference,  angular  cross-track,  and  angular  along-track  separations, 
respectively.  As  with  Figures  4-6,  the  first  plot  in  each  figure  shows  the  satellite  separations  described  by  the  truth 
orbits  and  the  COWPOKE  approximation,  and  the  second  plot  shows  the  differences  between  the  tmth  and 
COWPOKE  methods.  Thus,  the  second  plot  in  each  figure  is  the  error  inherent  in  the  given  COWPOKE  formulation 
for  this  highly  eccentric  test  case. 


Figure  10:  HEO  Test  Case  Radial  Differences  and  COWPOKE  Errors 


Figure  11:  HEO  Test  Case  Angular  Cross-Track  Differences  and  COWPOKE  Errors 
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Figure  12:  HEO  Test  Case  Angular  Along-Track  Differences  and  COWPOKE  Errors 


Figure  10  shows  the  COWPOKE  radial  error  to  be  significant  with  maximum  errors  being  on  the  same  order  of 
magnitude  as  the  relative  altitude  differences.  Figures  1 1  and  12  show  the  COWPOKJE  cross-track  and  along-track 
errors  are  on  the  order  of  10%.  The  sources  of  these  errors  can  be  determined  by  analyzing  the  approximations  used 
in  the  COWPOKE  formulation. 

Figure  13  plots  the  error  in  the  Fourier-Bessel  series  expansion  of  the  true  anomaly  of  the  reference  satellite.  From 
Figure  3,  we  expect  to  see  around  10%  error,  and  we  do  along  with  a  mixture  of  frequencies  in  the  signature  due  to 
the  series  truncation.  The  error  in  the  true  anomaly  approximation  maps  into  the  frequency  terms  and  becomes  the 
dominant  source  of  error  for  the  cross-track  component.  Figure  14  plots  the  error  in  the  approximation  of  the  true 
anomaly  difference;  the  error  in  this  approximation  is  an  order  of  magnitude  smaller  than  the  error  in  the  true 
anomaly  approximation  since  we  have  removed  an  order  of  eccentricity  in  the  differentiation  process  and  replaced  it 
by  an  eccentricity  difference.  Since  the  eccentricity  difference  is  an  order  of  magnitude  smaller  than  the  eccentricity, 
the  eccentricity  difference  error  is  much  smaller.  One  can  see  that  the  dominant  part  of  the  radial  and  along-track 
errors  is  due  to  the  COWPOKE  approximation  error  of  the  true  anomaly  difference  and  likely  the  eccentricity  series 
truncation  for  the  SM  frequency  terms.  If  one  chose,  these  errors  could  be  reduced  by  adding  higher  order  terms  to 
the  series  approximations  for  these  variables. 


•* 


Figure  13:  True  Anomaly  Approximation  Error  for  HEO  Test  Case 


Figure  14:  True  Anomaly  Difference  Approximation  Error  for  HEO  Test  Case 


Figure  15  plots  the  errors  in  the  COWPOKE  coordinate  system  approximations  for  the  HEO  test  case.  There  is  no 
plot  for  the  radial  conqionent  since  the  radial  approximation  was  not  used,  and  thus  the  coordinate  system  error  is 
zero.  The  cross-track  and  along-track  components  remain  around  or  below  the  1%  level  even  for  the  high 
eccentricity  case.  It  would  seem  that  the  coordinate  system  gives  the  COWPOKE  equations  an  error  bound  of  about 
1%. 
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Figure  15:  COWPOKE  Coordinate  System  Approximation  Errors  for  LEO  Test  Case 


The  cross-track  and  along-track  errors  for  the  HEO  case  are  in  line  with  expectations  given  the  series  truncations 
made  in  the  COWPOKE  formulation  for  this  test  case.  The  source  of  the  radial  error  was  larger  than  expected  or 
desired.  Further  analysis  of  Eq.  (4)  shows  that  there  is  coupling  between  the  true  anomaly  error  and  the  eccentricity 
difference.  The  only  way  to  mitigate  this  error  is  to  reduce  the  error  in  the  true  anomaly  difference  approximation 
which  can  be  accomplished  by  adding  additional  terms  to  the  series  expansion  or  by  solving  Kepler’s  equation 
iteratively.  This  would  also  serve  to  reduce  the  along-track  errors  as  well. 

Despite  the  errors  present  in  the  HEO  test  case,  the  COWPOKE  equations  still  model  bulk  of  the  relative  motion 
fairly  well,  and  the  geometric  foimdation  on  which  the  equations  are  built  can  provide  a  fair  amoimt  of  physical 
insight  into  the  relative  motion.  For  the  HEO  case,  the  physical  insight  is  very  similar  to  the  LEO  test  case  except  we 
now  observe  that  the  true  anomaly  does  not  vary  linearly  with  time.  While  the  series  expansions  do  not  provide  a 
great  deal  of  insight  into  how  the  true  anomaly  varies,  those  who  are  familiar  with  the  concept  should  be  able  to 
anticipate  the  trends.  In  Figures  10-12,  we  see  that  as  the  satellite  moves  towards  apogee  (recall  the  initial  mean 
anomaly  was  90  deg),  the  true  anomaly  rate  slows  which  pushes  the  relative  motion  extrema  to  toward  the  time  of 
perigee  passage.  It  is  interesting  to  note  that  initial  true  anomaly  difference  is  a  function  of  the  mean  anomaly 
difference,  eccentricity,  and  the  initial  mean  anomaly;  a  small  mean  anomaly  difference  at  perigee  of  a  highly 
eccentric  orbit  maps  into  a  much  larger  true  anomaly  difference  than  if  the  initial  conjjitions  occurred  at  apogee. 
Additionally,  one  must  be  aware  of  this  fact  when  considering  error  sources  in  the  COWPOKE  formulation  since 
only  first  order  differences  in  the  Keplerian  elements  have  been  included  in  the  equations  above.  Finally,  from  the 
COWPOKE  equations,  one  can  see  that  the  argument  of  perigee  also  plays  a  significant  role  in  the  relative  motion 


since  the  cross-track  motion  is  dependent  on  bofli  the  true  anomaly  and  the  argument  of  perigee.  The  argument  of 
perigee  essentially  controls  the  phasing  between  the  cross-track  and  radial/along-track  conqjonents  of  motion. 

CONCLUSIONS  AND  FUTURE  WORK 

This  paper  derived  the  Cluster  Orbits  With  Perturbations  Of  Keplerian  Elements  (COWPOKE)  equations  for 
xmperturbed  satellite  motion.  A  general  framework  is  provided  to  generate  a  set  of  equations  which  describe  the 
relative  motion  between  satellites  in  eccentric  orbits  explicitly  as  a  function  of  the  Keplerian  elements  of  the 
reference  satellite,  Keplerian  element  differences,  and  time.  This  is  acconqjlished  using  a  geometric  description  of 
the  separation  between  the  satellites  to  form  the  basis  of  motion  and  then  by  using  a  Fourier-Bessel  series  expansion 
of  the  true  anomaly  in  terms  of  the  mean  anomaly.  True  anomaly  and  radial  differences  are  derived  using  first  order 
perturbation  methods.  The  resulting  equations  are  meant  to  provide  accurate  representations  of  the  relative  motion 
between  satellites,  to  be  simple  to  inclement,  and  to  provide  physical  insight  into  the  relative  motion.  Test  cases 
show  that  this  has  been  accoiiq)lished  to  a  certain  degree. 

The  geometric  basis  on  which  the  COWPOKE  equations  are  built  is  accurate  to  the  1%  level  when  compared  to  the 
separation  distances.  This  is  currently  the  theoretical  accuracy  limitation.  For  most  applications  where  the 
eccentricities  are  below  0.1,  the  COWPOKE  equations  are  easily  realizable  and  can  be  as  accurate  as  the  hmitations 
of  the  geometric  basis;  however,  the  practical  limitation  for  highly  eccentric  cases  comes  from  generating  enough 
terms  in  the  series  expansion  of  the  true  anomaly  to  reduce  approximation  errors  to  an  acceptable  level.  This  work 
shows  that  the  radial  conq)onent  of  the  relative  motion  in  our  spherical  reference  system  is  particularly  sensitive  to 
true  anomaly  errors  for  high  eccentricity  orbits. 

The  next  step  in  the  development  of  the  COWPOKE  equations  will  be  to  include  dynamic  perturbation  effects  into 
the  relative  motion.  Reference  8  has  shown  how  secular  effects  due  to  Earth’s  oblateness  can  be  incorporated  using 
Keplerian  elements  and  element  differences.  Simulations  will  also  be  performed  to  determine  what  other 
perturbations  are  required  to  support  applications  such  as  modeling  geosynchronous  clusters.  We  will  also  look  into 
die  invertibility  of  flie  COWPOKE  equations.  Here  we  will  attempt  to  solve  for  Keplerian  element  differences  based 
on  desired  or  observed  relative  motion.  This  would  be  tremendously  useful  for  formation  design  and 
geosynchronous  close  approach  calculations. 
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ABSTRACT 

Optical  satellite  tracking  often  reveals  multiple  satellites  in  a  single  telescope  field  of  view.  The  Cluster  Orbits  With 
Perturbations  Of  Keplerian  Elements  (COWPOKE)  equations  are  used  to  estimate  the  relative  motion  of 
geosynchronous  satellites  and  determine  if  the  satellites  can  be  later  identified  based  solely  on  relative  position.  This 
paper  provides  the  development  of  the  COWPOKE  equations  for  modeling  the  relative  motion  of  geosynchronous 
satellites  and  anal3reis  demonstrating  the  feasibility  of  using  these  techniques  for  object  correlation.  Real  data 
relative  orbit  determination  results  are  provided  using  the  optical  tracking  assets  of  the  Air  Force  Maui  Optical  and 
Superconqiuting  (AMOS)  site. 


1.  INTRODUCTION 

Clusters  of  spacecraft  in  geosynchronous  orbit  (GEO)  are  becoming  more  common.  Orbital  slot  allocations  in  GEO 
are  rapidly  being  filled,  and  it  is  increasingly  difficult  to  acquire  slots  for  new  satellites.  Consequently,  many 
organizations  opt  to  collocate  flieir  spacecraft  in  the  same  slot.  Eutelsat  had  as  many  as  five  satellites  collocated  in  a 
formation  at  13°  E  by  2001  ['].  In  addition,  unintentional  close  approaches  have  occurred  and  could  become  more 
frequent  as  the  GEO  belt  becomes  more  crowded.  In  1997,  Telstar  401,  a  satellite  in  GEO,  experienced  a  major 
failure,  and  control  of  the  spacecraft  was  lost  [“].  It  is  drifting  in  GEO  and  has  had  encounters  with  other  satellites  as 
close  as  4  km  [‘“].  Since  that  time,  other  GEO  spacecraft  have  also  begun  drifting  uncontrollably.  Close  approaches 
and  satellite  formations  create  a  challenge  for  space  sinveillance.  Identification  of  the  satellites  in  these  clusters  can 
be  difficult,  and  cross-tagging  (misidentification)  occurs  n.  Resolvable  imaging  can  be  used  to  identify  spacecraft 
in  low  Earth  orbit  but  is  not  currently  possible  for  GEO  altitudes.  Non-imaging  approaches  include  comparing  the 
brightness  and  characterizing  the  color  of  light  reflected  from  the  spacecraft  for  identification  purposes  f]. 
However,  identifying  spacecraft  solely  from  their  dynamics  would  eluninate  flie  need  for  special  filters  or  other 
sensors.  A  better  understanding  of  the  relative  motion  of  the  spacecraft  could  reduce  cross-tagging  and  inqjrove 
close  approach  predictions.  In:5)roved  determination  of  minimum  approach  distances  could  eliminate  unnecessary 
collision  avoidance  maneuvers  and  minimize  propellant  usage. 

Cluster  orbits  can  be  modeled  using  the  Cluster  Orbits  With  Perturbations  Of  Keplerian  Elements  (COWPOKE) 
equations  of  motion  ["].  The  COWPOKE  equations  predict  spacecraft  separations  in  the  spherical  radial,  along- 
track,  and  cross-track  coordinate  frame  based  on  the  Keplerian  elements  of  the  reference  satellite,  the  element 
differences  for  the  second  satellite,  and  elapsed  time.  These  inputs  for  COWPOKE  can  be  obtained  with  relative 
metric  data  from  optical  sensors  or  space  surveillance  products.  Charge-coupled  device  (CCD)  images  of  clusters  of 
spacecraft  should  provide  very  accurate  measurements  of  spacecraft  separation  since  in-frame  error  sources 
theoretically  cancel  out.  The  United  States  Air  Force  currently  operates  several  optical  systems  capable  of  imaging 
clusters  at  the  Maui  Space  Surveillance  Complex,  such  as  Raven,  Phoenix,  and  the  Groimd-Based  Electro-Optical 
Deep  Space  Surveillance  System  (GEODSS).  Raven  uses  small,  commercially  available  telescopes  to  acquire  CCD 
images  of  space  objects  for  tracking  [™].  Astrometry  is  used  to  match  stars  in  the  CCD  image  to  the  star  catalog. 
Thus,  the  pointing  accuracy  of  the  images  can  approach  the  acciuacy  of  the  star  catalog  used.  Phoenix  is  a  Baker- 
Nunn  telescope  refurbished  to  take  wide-angle  CCD  images  [™‘].  As  many  as  21  satellites  have  been  detected  in  one 
Phoenix  image.  After  the  Deep  Stare  upgrade,  GEODSS  could  be  used  to  acquire  relative  metrics  as  well. 


This  paper  is  an  extension  of  previous  work  exploring  the  application  of  COWPOKE  towards  geosynchronous 
cluster  orbit  prediction  [“].  First,  a  perturbation  study  is  conducted  to  detemiine  the  force  modeling  required  at 
GEO.  Improvements  are  made  to  the  COWPOKE  equations  that  allow  for  better  representation  of  GEO  motion.  A 
method  is  constraed  to  estimate  the  Keplerian  element  differences  using  optical  measurements  of  relative  right 
ascension  and  declination.  Finally,  results  are  discussed  for  which  COWPOKE  was  used  to  predict  die  relative 
motion  of  a  cluster  of  satellites  in  GEO. 


2.  USING  COWPOKE  AT  GEO 

A  study  was  conducted  to  identify  perturbing  forces  that  significantly  affect  the  relative  motion  of  a  cluster  of 
geosynchronous  satellites.  The  effects  of  central  body  gravity,  third  body  gravity,  and  solar  radiation  pressure  were 
investigated  using  the  Draper  Semianalytic  Satellite  Theory  (DSST).  The  Draper  research  and  development  version 
of  the  Goddard  Trajectory  Determination  System  (DGTDS)  was  used  to  propagate  the  orbits  with  DSST  [*].  The 
motion  of  two  spacecraft  was  propagated  to  generate  a  “truth”  relative  motion  using  an  8x8  JGM-2  geopotential 
field,  luni-solar  third-body  point-mass  using  JPL  ephemerides,  and  solar  radiation  pressure  (SRP)  based  on  a 
spherical  satellite  and  cylindrical  Earth-shadow  model.  Next,  the  relative  motion  of  the  two  spacecraft  was 
generated  while  neglecting  one  of  the  sources  of  perturbations,  such  as  third-body  gravity.  The  radial,  along-track, 
and  cross-track  separations  obtained  from  the  “trath”  relative  motion  were  compared  to  those  generated  with  the 
inconplete  force  model  to  find  the  approximate  error  that  would  result  from  neglecting  that  perturbation  source  in 
the  COWPOKE  formulation.  These  test  runs  were  conducted  two  times;  the  first  set  of  runs  included  long  period 
effects  propagated  over  30  days,  and  the  second  set  investigated  short  periodic  effects  propagated  over  5  days. 
Table  1  shows  the  reference  elements  and  the  differential  elements  used  in  the  propagation. 


Table  1.  Initial  conditions  for  DSST  runs 
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Cases  that  neglected  all  non-spherical  gravity  forces  were  conducted  at  differing  initial  mean  anomalies.  Those 
omitting  highpr  order  geopotential  terms,  but  including  J2,  were  conducted  with  varying  longitudes  of  the  node  to 
siuvey  the  longitudinal  dependencies  of  the  tesseral  harmonics.  Tests  examining  luni-solar  effects  ran  at  varying 
days  of  fixe  year  while  those  for  SRP  effects  used  different  days  of  the  year  as  well  as  at  various  differential  area-to- 
mass  ratios  for  the  spacecraft.  Errors  were  calculated  by  dividing  the  difference  in  radial,  along-track,  or  cross-track 
separation  by  the  maximum  separations.  Fig.  1  shows  the  worst  case  errors  produced  by  neglecting  each  of  the 
perturbing  forces  after  30  days. 

The  SRP  results  in  Fig.  1  arise  from  a  case  in  which  the  area-to-mass  ratio  of  one  satellite  is  ten  times  that  of  the 
other.  As  the  separation  between  the  spacecraft  decreases  to  zero,  fire  effect  of  neglecting  gravitational  effects  also 
decreases  to  zero;  however,  SRP  effects  do  not  decrease  to  zero  if  the  area-to-mass  ratios  differ.  Because  of  this, 
SRP  effects  should  be  taken  into  account  if  high  accuracy  is  warranted  over  long  periods  or  if  the  satellites  are  within 
several  kilometers. 

The  effects  of  SRP  on  an  orbit  were  formulated  using  Keplerian  elements  and  the  element  differences,  assuming  that 
the  latter  quantities  are  small.  However,  because  Keplerian  elements  are  singular  for  i  =  0  and  e  =  0,  orbits  in  GEO 
can  have  large  differences  in  Q,  o),  and  Mo  and  still  remain  wdfiiin  a  few  kilometers  of  each  other.  Simulations  found 
that  this  formulation  of  SRP  did  not  prove  to  be  useful  for  GEO  and  therefore  was  not  used  in  later  tests.  However, 
for  most  cases,  two-body  motion  results  in  acceptably  low  error. 
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Fig.  1.  Maximum  error  encountered  in  the  separation  of  two  satellites  at  GEO  due  to  the  neglect  of  a 
perturbing  force  after  a  30-day  orbit  propagation. 

Another  source  of  orbit  perturbation  is  stationkeeping  and  momentum  control  maneuvers.  Since  maneuvers  were  not 
modeled  in  COWPOKE,  they  present  an  additional  source  of  error  if  occurring  during  the  observation  period. 

3.  IMPROVED  COWPOKE  EQUATIONS 


Alfliough  COWPOKE  has  been  shown  to  be  an  effective  predictor  of  the  relative  motion  of  a  cluster  of  satellites, 
there  are  some  sources  of  error.  For  GEO,  the  right  ascension  of  the  ascending  node,  argument  of  perigee,  and  mean 
anomaly  element  differences  may  not  be  small  which  violates  the  assumptions  of  the  original  COWPOKE  derivation. 
In  particular,  the  approximation  of  the  true  anomaly  difference  was  linear  in  terms  of  the  mean  anomaly  difference; 
this  causes  error,  particularly  in  the  along-track  direction,  when  the  mean  anomaly  difference  is  significant.  The 
cross-track  term  conqronent  of  COWPOKE  also  showed  error  with  large  right  ascension  differences. 


Replacing  the  linear  approximation  of  the  true  anomaly  difference  with  an  exact  difference  of  the  two  true  anomaly 
terms  substantially  reduced  the  along-track  error.  For  the  near-circular  GEO  case,  the  true  anomaly  of  each  satellite 
as  a  function  of  time  is  approximated  by  a  first  order  expansion  in  terms  of  eccentricity  and  mean  anomaly  which 
was  shown  to  be  effective  in  [6].  The  cross-track  corrqjonent  was  improved  by  an  investigation  into  the  spherical 
geometry  involved.  Corrections  were  made  to  both  amplitude  and  phasing  of  the  cross-track  corr^onent  which  are 
accurate  for  large  values  of  right  ascension  of  the  ascending  node,  argument  of  perigee,  and  mean  anorrraly 
differences.  The  irtqrroved  COWPOKE  equations  are 
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Sv  =  SM  +  2(e  +  <5e)sin(M  +  SM ) — 2esin(M ) 


Sr  is  the  separation  in  the  radial  direction,  Sxi  is  the  cross-track  separation,  and  is  the  along-track  separation. 
5v  is  the  difference  in  true  anomaly,  a,  e,  i,  Q,  a>,  and  M  are  the  orbital  elements  of  the  reference  satellite,  and  Sa , 
Se,  Si,  S£2,  5m,  and  ^are  the  differences  in  the  elements  of  the  two  satellites.  JU  is  the  gravitational  parameter,  and 

t  is  the  time  elapsed  since  the  epoch  of  M . 


The  along-track  term  still  exhibits  error  when  the  8Q  and  5i  terms  are  large.  To  avoid  this  as  an  error  soince  in  the 
relative  orbit  determination  experiments,  it  was  decided  to  use  the  along-track  component  from  Vadali’s  umt  sphere 
model  for  relative  motion  [*'].  The  Vadali  imit  sphere  model  is  very  similar  in  philosophy  to  COWPOKE;  relative 
motion  is  modeled  through  Keplerian  elements  and  element  differences.  Vadali’s  geometric  model,  however,  is 
more  rigorous  than  the  single  COWPOKE  approach.  Here  is  the  along-track  component  of  Vadali’s  unit  sphere 
model  which  was  incorporated  into  COWPOKE  for  the  remainder  of  this  analysis: 


—  =  cos^(f  /2)  cos^  ((/  -b  Si)  IT)  sin(^ry  +  Sv  +  SOi) 
r 

-1-  sin^  (i  /2)  sin^  ((i  -i-  Si)l2)  sin(^0  +  Sy-  SCi) 

—  sm(i/2)cos^((i+  Si)  /2)sin{2{co  +  v)  +  S(o  +  Sy+ 

-  cos^(/  /2)sin^((i  -I-  Si)  /2)cos(2((»  +  y)  + So  +  Sy-S£i£) 


+  ^  sin(z)  sin(/  +  [sin(<5to  +  Sy)  +  sin(2(<»  +  y)-\-5o}-\-  ^v)] 


(3) 


4.  DETERMINATION  OF  ELEMENT  DIFFERENCES 

Predicting  the  relative  motion  of  a  cluster  of  satelhtes  with  COWPOKE  requires  a  set  of  orbital  elements  for  the 
reference  satellite  and  the  relative  elements  of  the  second  satellite.  A  least-squares  orbit  determination  method  was 
used  to  find  these  element  differences  [*"]. 

COWPOKE  expresses  spacecraft  separations  in  the  spherical  radial,  along-track  and  cross-track  reference  frame,  but 
the  optical  observations  used  in  this  effort  are  in  the  topocentric  right  ascension  and  declination  frame.  Therefore, 
the  topocentric  observations  had  to  be  converted  to  geocentric  observations,  which  requires  satellite  range 
knowledge  as  well  as  the  local  sidereal  time  For  this  analysis,  a  constant  range  value  was  used  and  values  of 
UTl-UTC,  precession  and  nutation  angles,  and  lunar  terms  were  ignored;  it  is  believed  that  these  approximations  do 
not  have  significant  impact  once  the  observations  are  differenced.  The  COWPOKE  cross-track  and  along-track 
separations  could  then  be  equated  with  the  geocentric  right  ascension  and  declination  frame  as  follows: 

5a  =  SatcosO  -  sin^ 

(4) 

Sd  =  Sxt  COS0  4-  Sat  sin  0 


where 

0  =  fcos(fi>-l-u)  (5) 

Let  y  be  the  relative  observation  vector,  and  X  the  state  vector  containing  the  orbital  element  differences.  Eqs.  (4)  & 
(5)  represent  a  nonlinear  mapping  between  die  state  vector  and  the  observations.  In  order  to  invert  the  problem,  we 
must  linearize  the  equations  about  a  reference  trajectory,  X*. 


(6) 


7  =  M  =  F(X)«F(X*)  +  —  (X-X*) 

_Sd\  aAT;,. 

X  =  lSa  Se  Si  tSQ  da  SMJ  (7) 

The  ^ateims  are  the  observed  differences  in  right  ascension,  and  the  5d  are  the  observed  differences  in  declination. 
Eq.  (6)  can  be  rearranged  and  terms  can  be  redefined  as  follows: 

(x-z*) 

y^Hx  (8) 

where  y  is  the  difference  between  the  observed  and  calculated  relative  observations,  H  is  the  linearized  observation- 
state  relationship,  and  x  is  the  estimated  correction  to  the  state  matrix.  If  at  least  3  relative  observation  pairs  are 
included  in  the  jv  vector,  one  can  estimate  the  state  deviation  as  shown  below: 

x  =  {H^HyH^y  (9) 

The  estimate  of  the  state  can  be  updated  in  an  iterative  fashion,  as  shown  below,  imtil  the  solution  converges. 

xl,=x':+x  (10) 

Using  this  method,  an  estimate  of  the  Keplerian  element  differences  was  obtained.  This  method  requires  tiiat  the 

Keplerian  elements  of  the  reference  satellite  are  known.  The  NORAD  two-line  element  set  (TLE)  of  the  reference 
satellite  was  used  for  the  reference  orbit  in  this  study.  With  the  element  differences,  the  relative  right  ascension  and 
declination  of  the  two  satellites  can  be  predicted  using  COWPOKE. 

5.  SIMULATION  STUDY 

In  order  to  test  the  feasibility  of  using  COWPOKE  to  better  predict  relative  motion,  a  simulation  was  performed 
using  2  collocated  geos5mchronous  satellites.  Truth  orbits  were  propagated  using  the  Cowell  Special  Perturbations 
(SP)  propagator  internal  to  DGTDS.  The  truth  orbits  spanned  from  Jan  1  -  Feb  5,  2003.  Table  2  contains  the 
ngnilating  orbital  elements  for  the  satellites  referenced  to  the  mean  equator  and  mean  equinox  of  the  B  1950.0 
coordinate  system. 

Table  2.  The  initial  orbital  elements  of  the  two  satellites  used  in  the  identification  study 
_ with  epoch  January  1, 2003, 0  hr. _ 


Satellite  1 

Satellite2 

Semimajor  Axis  (km) 

42164 

42164 

Eccentricity 

0.000512 

0.000812 

Inclination  (deg) 

0.05 

0.05 

Right  Ascension  of  the 

Ascending  Node  (deg) 

20 

340 

Argument  of  Perigee  (deg) 

33 

67 

Mean  Anomaly  (deg) 

0 

6 

The  next  step  in  the  process  was  to  develop  TLE  representations  of  the  truth  orbits.  This  step  was  necessary  since 
TLEs  are  used  in  the  observation  correlation  process.  To  do  this,  an  orbit  using  the  GP4/DP4  propagator  internal  to 
DGTDS  was  fit  to  the  position  and  velocity  vectors  produced  by  the  Cowell  ttuth  trajectory.  The  position  and 


velocity  data  spanned  1  Jan  to  1  Feb  and  were  spaced  at  every  hour.  The  resulting  fits  were  acciuate  to  around  the  2 
km  level  RMS.  Fig.  2  plots  the  DP4  trajectories  relative  to  the  troth  orbit  of  Sat  1  over  February  1-5,  a  5-day 
prediction  interval. 

One  can  clearly  see  in  Fig.  2  that  the  TLE  trajectory  for  Sat  2  comes  closer  to  the  Troth  location  for  Sat  1  than  the 
Troth  location  for  Sat  2  during  a  significant  portion  of  the  orbit.  Similarly,  the  TLE  trajectory  for  Sat  1  comes  close 
to  the  Troth  location  for  Sat  2  during  a  portion  of  the  orbit.  Both  of  these  situations  could  lead  to  a  cross-tag  in  the 
observation  correlation  process  (i.e.,  observations  of  Sat  2  could  be  incorrectly  attributed  to  Sat  1  and  vice  versa). 


Fig.  2.  TLE  orbit  predictions  of  the  motion  of  Satellite  2  relative  to  Satellite  1  truth  orbit 
compared  to  the  truth  relative  motion. 

The  tnffan  orbital  elements  firom  the  TLE  fits  (epoch  1  Feb)  were  differenced  and  used  to  initialize  the  COWPOKE 
equations  along  with  the  TLE  elements  for  Sat  1.  One  can  see  in  Fig.  3  that  even  using  the  flawed  initial  conditions 
produced  by  the  TLE  fits,  the  relative  motion  is  very  representative  of  the  Troth  relative  motion.  This  signifies  that 
the  relative  position  is  a  powerful  piece  of  information  that  can  be  used  to  help  eorrelate  optical  observations  at  the 
sensor. 


Fig.  3.  COWPOKE  relative  motion  predictions  with  TLE  and  truth  orbits. 


6.  REAL  DATA  RESULTS 

Observations  from  Raven  were  used  to  test  COWPOKE’S  effectiveness.  Raven  images  were  taken  of  die  DirecTV 
4S  and  AMC  4  spacecraft  collocated  at  101°  West  longitude  during  the  nights  of  July  23-24  and  July  29-August  1, 
2003.  One  of  these  images  is  shown  in  Fig.  4. 

Raven  images  were  used  to  con:5)ute  the  separation  of  the  two  satellites.  Sat  1,  the  reference  satellite,  was  chosen  to 
be  AMC-4,  and  DirecTV  4S  was  designated  Sat  2.  The  reference  orbit  for  Sat  1  was  generated  using  the  TLE  from 
July  20,  2003.  Other  TLEs  might  be  available  at  an  epoch  closer  to  the  24*  but  for  R&D  Raven  operations,  the 
catalog  is  only  updated  every  few  days.  The  observed  separation  in  right  ascension  and  declination  on  the  night  of 
the  23'**  were  used  to  estimate  the  Keplerian  element  differences.  The  observations  from  the  23''*  didn’t  span  a  long 
enough  period  to  accurately  solve  for  the  difference  in  semimajor  axis,  so  an  a  priori  estimate  of  0  m  was  added  for 
Sa,  with  a  standard  deviation  of  1000  km.  With  the  resulting  estimate  of  the  element  differences,  COWPOKE  was 
used  to  predict  the  relative  position  of  the  two  satellites  at  the  time  of  each  observation  taken  on  the  night  of  the  24*. 
The  COWPOKE  predictions  are  con^ared  to  the  positions  predicted  by  the  TLEs  at  flie  time  of  each  observation  in 
Fig.  5.  The  COWPOKE  prediction  of  the  location  of  Sat  2  was  off  by  an  average  of  155  microradians,  while  the 
TLE  predictions  differed  from  the  observations  by  an  average  of  576  microradians. 


Fig.  4.  Raven  image  of  AMC-4  and  DirecTV  4S 


Motion  Relative  to  Sat  1  Truth  Orbit  on  7-24-03 
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Fig.  5.  COWPOKE’S  prediction  of  the  position  of  Satellite  2  on  July  24  compared 
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The  observations  from  July  23  and  24  were  used  to  estimate  the  element  differences  and  predict  the  relative  motion 
for  July  29,  when  the  next  telescope  images  were  taken.  DP4  predictions  and  the  COWPOKE  reference  orbit  were 
obtained  using  TLEs  for  Sat  1  and  Sat  2  from  the  26“"  and  27*  respectively.  The  results  for  the  29*  are  shown  in 
Fig.  6.  The  COWPOKE  predictions  differed  from  the  truth  by  an  average  of  390  microradians,  while  the  TLE 
predictions  were  off  by  721  microradians. 


Motion  Relative  to  Sat  1  Truth  Orbit  on  7-29-03 
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Fig.  6.  COWPOKE’S  prediction  of  the  position  of  DirecTV  4S  on  July  29  compared  to  the  TLE  predictions 

The  observations  from  July  23, 24,  and  29  were  then  used  to  estimate  improved  element  differences.  Those  element 
differences  were  used  to  predict  the  relative  position  of  Sat  2  on  the  30*,  and  the  results  are  plotted  in  Fig.  7. 
COWPOKE  had  an  average  error  of  300  microradians,  and  the  TLE  for  Sat  2  averaged  894  microradians.  Fig.  8 
shows  the  results  of  similar  predictions  for  July  31  using  the  observations  from  all  previous  nights.  To  be  clear,  the 
TLE  prediction  span  at  this  point  is  several  days  while  the  COWPOKE  prediction  is  only  one  day. 


Fig.  7.  COWPOKE’S  prediction  of  the  position  of  DirecTV  4S  on  July  30  compared  to  the  TLE  predictions 

In  Fig.  8,  one  can  see  that  the  true  position  of  Sat  1  has  shifted  away  from  the  TLE  prediction  since  the  night  before. 
Also,  the  COWPOKE  prediction  no  longer  matches  the  observation  for  Sat  2,  especially  in  declination.  There  is 
strong  evidence  that  Sat  1  performed  a  stationkeeping  maneuver  between  the  30*  and  the  31““.  Even  with  the 


possible  maneuver,  COWPOKIE  still  provided  a  better  estimate  of  the  relative  motion  than  the  TLE.  COWPOKE 
had  an  average  error  of  450  microradians,  and  TLE  for  Sat  2  averaged  869  microradians. 


Whatever  the  cause  of  the  sudden  shift,  it  was  decided  to  start  a  new  fit  span.  Only  the  observations  from  the  night 
of  the  31“  were  used  to  solve  for  the  new  element  differences,  and  fiiese  were  propagated  forward  with  COWPOKE 
to  the  1“  of  August.  The  COWPOKE  and  TLE  predictions  for  that  night  are  compared  to  the  observations  in  Fig.  9. 
COWPOKE  had  an  average  error  of  210  microradians,  and  the  TLE  for  Sat  2  averaged  958  microradians.  One  of 
the  advantages  of  COWPOKE  is  that  the  effects  of  maneuvers  can  be  mitigated  by  using  a  one  day  fit  span  while 
TLEs  typically  use  a  longer  fit  span  and  will  suffer  the  effects  of  maneuvers  continuously. 


Fig.  9.  COWPOKE’S  prediction  of  the  position  of  DirecTV  4S  on  August  1  compared  to  the  TLE  predictions 


While  the  COWPOKE  results  are  somewhat  encouraging  when  coiiq)ared  to  the  TLE’s,  the  overall  performance  is 
not  as  good  as  expected.  Perturbation  analysis  indicate  that  the  equations  should  be  accurate  with  only  a  few  percent 
error.  The  simulation  results  showed  similar  error  levels.  If  the  relative  orbit  estimation  algorithms  were  functioning 
properly,  one  would  expect  to  see  results  with  errors  at  the  few  percent  level.  For  the  one  day  fit  cases  shown  in  Fig. 
5  and  9,  larger  errors  might  be  expected  due  to  limited  observability  over  a  short  data  arc.  For  fixe  five  day 
prediction  case  shown  in  Fig.  6,  one  might  also  expect  to  see  larger  errors  due  to  the  prediction  interval.  Then  there 
is  tile  case,  shown  in  Fig.  8,  where  the  reference  satellite  appears  to  maneuver;  this  would  also  cause  prediction  error. 
However,  one  case  remains,  shown  in  Figure  7,  where  several  days  of  data  are  used  in  the  estimation  process  and  the 
prediction  interval  is  only  one  orbit;  the  COWPOKE  prediction  error  is  as  large  as  all  of  tiie  other  cases  and  is 
around  10%  of  the  separation  distance.  This  is  larger  than  expected  and  indicates  that  there  may  be  an  unknown 
error  source  in  the  algorithm  or  software  tools.  Regardless,  efforts  must  be  made  to  better  understand  the  limitations 
of  this  approach. 


7.  CONCLUSIONS  AND  FUTURE  WORK 

This  work  has  shown  that  the  COWPOKE  equations  can  be  used  to  provide  meaningful  relative  motion  of 
geosynchronous  satellite  clusters.  Perturbation  analysis  indicated  that  2-body  dynamics  are  adequate  for  medium 
accuracy  applications.  Improvements  were  made  to  the  equations,  however,  to  account  for  large  right  ascension  of 
the  ascending  node,  argument  of  perigee,  and  trae  anomaly  element  differences. 

Estimating  the  Keplerian  element  differences  and  using  the  COWPOKE  equations  to  predict  the  relative  motion  can 
supply  valuable  information  in  spacecraft  identification.  Using  six  nights  of  Raven  images,  it  is  shown  that 
COWPOKE  estimated  the  position  of  DirecTV  4S  relative  to  AMC-4  much  better  than  TLE  predictions,  even  with 
unmodeled  maneuvers.  This  indicates  that  COWPOKE  holds  the  potential  to  be  a  valuable  space  surveillance  tool. 

Sizable  error  still  remains  in  the  relative  orbit  prediction  results.  These  errors  are  larger  than  anticipated  so  care 
must  be  taken  to  detemine  the  major  source  of  this  error  and  remove  it.  If  this  can  be  accon^lished,  the  relative  orbit 
estimation  approach  Avill  be  far  more  valuable.  Beyond  those  improvements,  flie  effects  of  differential  SRP  should 
be  formulated  for  circular,  equatorial  orbits  to  better  predict  the  motion  of  GEO  clusters. 
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The  motion  of  satellites  flying  in  formation  around  the  Earth  is  ordinarily  modeled  by 
HilPs  Equations  of  relative  motion.  Recently,  Sabol  et.  al.  developed  a  set  of  more  intuitive 
equations,  known  as  the  Cluster  Orbits  With  Perturbations  Of  Keplerian  Elements 
(COWPOKE)  Equations.  To  accurately  model  orbital  motion,  these  equations  must  be 
modifled  to  incorporate  various  perturbation  effects,  among  them  J2  and  drag.  These  effects 
have  been  added  to  the  model,  and  simulations  were  run  to  compare  this  approach  with 
"truth”  orbits  generated  by  Analytical  Graphics,  Inc.’s  Satellite  Tool  Kit. 
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-  denotes  differential 

=  radial  position,  km 

=  cross-track  position,  km 

=  along-track  position,  km 

=  semimajor  axis,  km 

=  eccentricity 

=  mean  anomaly,  rad 

=  true  anomaly,  rad 

=  cross-track  angular  separation 

=  along-track  angular  separation 

=  right  ascension  of  the  ascending  node,  rad 
=  argument  of  perigee,  rad 

=  inclination,  rad 

=  mean  motion,  1/sec 

==  semilatus  rectum,  km 

=  Earth  oblateness  parameter,  .0010826269 

Earth  radius,  6378.1363  km 
=  denotes  time  derivative 

=  Earth  gravitational  parameter,  398600.44 1 5  km^ /sec^ 

=  time,  seconds 

=  denotes  initial  condition 

=  denotes  reference  satellite 

=  denotes  second  satellite 

=  placeholder  for  mean  motion  derivative 

=  drag  force,  N 

=  velocity  magnitude,  m/sec 

=  ballistic  coefficient,  kg/m^ 

=  density,  kg/m^ 

=  mass,  kg 

=  drag  coefficient 


=  cross-sectional  area,  km^ 

=  scale  height,  km 

height  above  the  ellipsoid,  km 
==  reference  height,  km 

=  reference  density,  kg/m^ 


Introduction 


INCH  1878,  the  dominant  equations  of  motion  to  describe  relative  satellite  position  and  velocity  have  been  Hill's 
Equations.  Later  modified  by  Clohessy  and  Wiltshire,  these  equations  form  the  simplest,  most  basic  model  of 
satellite  formation  flight  in  use  today.  Of  course,  for  most  applications,  a  robust  orbit  propagator  such  as  the 
Goddard  Trajectory  Determination  System  (GTDS)  provides  highly  accurate,  near-truth  results.  However,  on  some 
occasions,  propagators  are  not  available,  or  an  analytical  representation  is  simply  preferred.  It  is  in  these  instances 
that  Hill's  Equations  are  ordinarily  used  to  model  formation  flight. 

Hill's  Equations,  in  their  simplest  form,  are  derived  directly  from  the  two-body  equations  of  relative  motion  and 
do  not  include  any  orbital  perturbations.  Although  it  is  possible  to  incorporate  perturbations,  the  equations  were 
intended  to  model  the  close  approach  of  two  satellites  over  a  short  period  of  time,  not  for  indefinite  periods  as  is  the 
case  for  formation  flight;  the  equations  tend  to  break  down  over  time.  Hill's  Equations  assume  that  the  orbit  is 
circular  and  that  the  satellites  are  close  together:  although  this  works  well  for  approach  and  rendezvous  scenarios,  the 
assurrq)tions  are  not  necessarily  valid  for  most  formation-flight  applications.  Also,  while  Hill's  Equations  follow 
naturally  from  the  algebraic  equations  of  motion  of  the  satellites,  the  approach  can  be  non-intuitive  at  best. 

Recently,  Sabol  et.  al.  described  a  novel  geometric  approach  to  formation  flight  modeling:  the  Cluster  Orbits 
With  Perturbations  Of  Keplerian  Elements,  or  COWPOKE,  Equations.  ‘  Rather  than  using  a  Cartesian,  satellite- 
centered  reference  frame,  the  COWPOKE  Equations  are  based  on  a  spherical  Earth-centered  frame  and  use  known 
Keplerian  orbital  element  differences  to  directly  determine  the  cross-track,  along-track,  and  radial  separations 
between  two  satellites.  Not  only  is  this  direct,  geometric,  Keplerian  approach  more  intuitive  than  the  algebraic  Hill 
model,  but  because  existing  perturbation  models  are  ordinarily  given  in  terms  of  Keplerian  element  differences, 
integrating  those  models  into  ^e  equations  of  motion  is  made  considerably  simpler. 

Prior  to  this  investigation,  however,  incorporation  of  any  perturbations  into  the  COWPOKE  Equations  had  not 
been  attempted;  although  Hill  et  al.  have  iir^roved  and  extended  the  equations  for  use  on  geosynchronous  orbit.^ 
Here,  I  describe  the  steps  taken  to  include  both  the  h  (Earth  oblateness)  and  atmospheric  drag  perturbations  into  the 
COWPOKE  Equations,  and  compare  the  resulting  propagations  with  those  generated  by  Satellite  Tool  Kit's  SGP4 
propagator.  I  have  also  included  in  the  Appendix  a  similar  derivation  for  Hill's  Equations  with  perturbations, 
although  these  results  have  not  been  correlated  with  the  COWPOKE  or  STK  propagations. 


Approach 

The  inproved  COWPOKE  Equations  are  as  follows: 
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where  dr  is  the  radial  displacement,  a  and  p  are  the  along-track  and  cross-track  angular  separations,  respectively, 
a  is  die  semimajor  axis  of  the  reference  satellite's  orbit,  e  is  its  eccentricity,  i  its  inclination,  co  its  argument  of 
perigee,  Q  its  right  ascension  of  the  ascending  node,  M  its  mean  anomaly,  v  its  true  anomaly,  and  r  die  radius.  5at 
and  5xt  denote  along-track  and  cross-track  linear  separations.  For  details  on  the  derivation  of  these  equations,  see 
Refs.  1  and  2. 

Addition  of  the  h  (first-order  Earth  oblateness)  perturbation  effects  secular  changes  in  the  right  ascension  of  the 
ascending  node,  argument  of  perigee,  and  mean  anomaly: 
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where  n  = 


is  the  mean  motion,  //  the  Earth's  gravitational  parameter,  p  the  semilatus  rectum,  J2  the  first- 


order  Earth  zonal  harmonic  coefficient,  and  the  Earth's  radius. 

An  inspection  of  Equations  1,  2,  and  3  will  reveal  that  inserting  these  secular  effects  into  the  COWPOKE 
Equations  requires  deriving  new  expressions  for  the  SQ,  5co^  and  8v  terms.  The  8Q  term  is  found  by  the  following 
method: 
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where  orbital  elements  of  satellite  2  are  expressed  as  the  elements  of  the  reference  satellite  plus  a  differential,  and 
the  subscript  o  denotes  an  initial  given  quantity.  Similarly, 
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The  true  anomaly  is  a  bit  more  difficult.  Trae  anomaly  difference  is  approximated  by  a  first  order  expansion  in 
terms  of  eccentricity  and  mean  anomaly:^ 

3v  =  3M  +  2{e  +  3e)  sin(M  +  3M)  -  2e  sin(M)  (5) 

Therefore,  to  determine  the  change  in  true  anomaly,  the  change  in  mean  anomaly  must  first  be  derived, 
no  perturbations,  M  =  fl ,  for  the  J2  case,  M  =  n  +  M ,  with  M from  Equation  2.  We  know  that  fl 
so  letting  CO^  =n  +  Mj^ ,  SM  can.  be  found  as  above: 
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All  of  these  new  element  differences  are  inserted  back  into  the  original  COWPOKE  Equations  (Eqn.  1),  and  if 
the  only  perturbation  desired  is  J2,  the  derivation  stops  here. 

Incorporation  of  the  atmospheric  drag  perturbation  follows  a  similar  course.  The  force  of  drag  is 
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with  p  the  atmospheric  density,  V  the  along-track  velocity,  and  BC,  the  ballistic  coefficient  of  die  satellite,  equal 


to 


where  tn  is  the  satellite's  mass,  Cd  is  its  drag  coefficient,  and  A  is  its  area  normal  to  the  velocity.  Atmospheric 
density  will,  of  course,  vary  throughout  the  satellite's  orbit.  There  are  a  number  of  approximations  in  use  to  calculate 
density;  I  have  selected  a  simple  exponential  model  for  purposes  of  this  simulation,  where 


^ellip 

p  =  Po~  n  (9) 

with  po  and  ho  reference  density  and  reference  altitude,  hoUip  the  actual  attitude  above  the  ellipsoid,  and  H  the 
scale  height.  These  quantities  can  be  found  in  a  table  such  as  Table  8-4  in  Reference  3. 

Given  Equations  7,  8,  and  9,  the  variation  of  parameters  method  is  used  to  determine  the  element  differences  due 
to  drag:^ 
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Again,  as  above,  die  final  Keplerian  element  differences  become 
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where,  again,  the  subscript  o  denotes  an  original  given  element  difference,  and  dcoj2  is  the  argument  of  perigee 
difference  due  to  Earth  oblateness  as  found  above.  These  final  element  differences  are  substituted  back  into 
Equations  1,  and  the  derivation  of  the  COWPOKE  Equations  with  J2  and  atmospheric  drag  perturbations  is  con5)lete. 
An  examination  of  the  Appendix  should  convince  the  reader  that  this  method  is  considerably  simpler  than  the 
incorporation  of  perturbations  into  Hill's  Equations. 


Methodology  and  Results 

To  validate  and  test  the  above  derivations,  simulations  were  run  using  a  set  of  mean  orbital  elements  with  all 
element  differences  equal  to  .01.  The  elements  and  differences  are  given  in  Table  1. 

The  modified  COWPOKE  Equations  are  modeled  in  Matlab  as  above.  The  simulations  run  for  60  hours  with 
data  points  taken  at  1 -minute  intervals.  To  validate  and  verify  the  results,  simulations  are  also  run  using  STK's  SGP4 
propagator.  In  STK,  the  relative  position  of  the  two  satellites  is  automatically  calculated  as  a  vector;  the  data  is 
exported  to  Matlab  and  compared  with  the  COWPOKE-derived  orbit.  For  the  drag  perturbation,  the  mass  of  both 
satellites  is  taken  to  be  1000  kg  and  the  cross-sectional  area  .000001  km^.  The  drag  coefficient  of  satellite  1  is  2,  and 
the  drag  coefficient  of  the  second  satellite  is  slightly  varied,  at  2.001.  This  results  in  an  approximate  250000  kg/km^ 
difference  in  ballistic  coefficient.  These  parameters  are  shown  in  Table  2. 

Figures  1  through  3  show  the  radial,  along-track,  and  cross-track  differences,  respectively,  for  both  propagation 
methods;  the  second  plot  in  each  figure  shows  the  error  between  COWPOKE  and  STK. 


Table  1:  Keplerian  Elements  and  Element  Differences 

Reference  Element 

_ Elements _ Differences 

a  (km)  T700Q _ .01 

e  .01  .01 


Sat  1  drag  2 


coefficient _ 

Sat  2  drag  2.001 

coefficient _ 


Conclusion 

The  COWPOKE  Equations  have  been  modified  to  include  the  effects  of  both  zonal  harmonic  and  atmospheric 
drag  perturbations.  This  was  expected  to  increase  the  accuracy  of  the  equations  wifli  respect  to  a  "truth"  orbit  as 
propagated  by  a  commercial  orbit  integration  program,  rendering  them  more  useful  as  a  simple  model  of  relative 
motion  between  satellites  in  formation  flight  scenarios.  When  compared  to  die  orbit  produced  by  STK's  SGP4 
propagator,  the  modified  COWPOKE  equations  instead  showed  somewhat  decreased  accuracy.  A  number  of  reasons 
can  be  put  forth  for  this,  notably  the  different  propagation  methods  used  in  this  and  previous  papers.  Before  drawing 
a  final  conclusion  about  the  modified  COWPOKE  Equations  as  here  derived,  a  similar  simulation  should  be  run 
utilising  the  GTDS  software  -  the  same  propagator  t^t  was  used  in  the  original  COWPOKE  papers.^  It  would 
similarly  be  constractive  to  conqiare  both  the  COWPOKE  and  STK  results  with  file  orbit  defined  by  Hill's  Equations 
with  perturbations,  which  are  derived  in  the  Appendix.  Furthermore,  simulations  should  be  run  using  observational 
data  from  real  satellites  to  conqiare  the  COWPOKE  results  wifti  reality.  If,  indeed,  such  research  as  indicated  above 
should  prove  these  modified  COWPOKE  Equations  accurate  and  reliable,  additional  perturbations  such  as  solar 
radiation  pressure  and  third-body  gravitational  effects  can  be  incorporated  into  the  equations,  further  constraining  the 
forces  acting  on  the  satellite  and  drawing  the  simulation  ever  closer  to  the  actual,  "tmth"  position  of  the  orbit. 

Appendix:  Hill's  Equations  Solved  with  Ji  and  Drag  Perturbations 

Hill's  Equations  solve  the  relative  two-body  equation  of  motion: 
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for  two  satellites  traveling  close  to  each  other  on  a  near-circular  orbit.  The  equations  are  as  follows: 
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where  the^  are  components  of  .  The  equations  are  solved  many  places,  including  Ref.  3;  the  position  and 
velocity  solutions  are: 
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The  secular  effects  of  &st-order  zonal  harmonics  (J2)  on  co,  Q,  and  Mare  given  in  Equations  2.  McLaughlin  et. 
al.*  have  shown  die  method  for  introducing  these  secular  variations  into  Hill's  Equations,  along  with  differential 
effects  due  to  inclination  differences  between  the  two  satellites. 

Adding  the  drag  perturbation  to  Hill's  Equations  is  more  complicated.  In  this  drag-and-J2-only  case,^  and^  in 
Equation  13  are  both  zero,  and^  is  equal  to  Fo,  since  the  drag  force  acts  only  in  opposition  to  the  along-track 
velocity  vector.  Starting  with  Equation  7,  a  relative  drag  force  is  calculated  by  differencing  the  forces  for  the  two 
satellites.  Each  individual  parameter  for  the  second  satellite  is  expressed  as  the  value  for  satellite  1  plus  a 
differential.  For  example,  V2  is  expressed  as 


Fj  =V,+SV 
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.■Simplifying  and  removing  second-order  and  higher  differentials  from  the  resulting  expression,  the  relative  drag 
force  can  finally  be  expressed  as 
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Setting^  in  Equation  13  equal  to  fD,  and  recognizing  that  the  z  force  is  uncoupled  to  the  x  andy  conq)onents,  the 
system  of  two  equations  can  be  solved.  The  solution  process  is  as  follows:  Differentiate  the  x  equation  once,  then 
solve  via  Laplace  transforms,  and  differentiate  the  result  to  arrive  and  the  x-velocity  solution.  The  y  equation  can  be 
solved  by  simple  integrate.  The  resulting  drag  solution  is  then  merged  with  the  J2  solution^  to  arrive  at  the  following 
set  of  equations  describing  the  motion  of  one  satellite  with  respect  to  another: 
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Despite  the  rather  perfunctory  treatment  given  die  derivation  in  this  Appendix,  it  should  be  clear  from  a  glance  at 
the  preceding  equations  that  working  with  the  simple  elegance  of  the  COWPOKE  Equations  is  eminently  preferable 
to  tWs  mess,  in  addition  to  providing  an  inherently  more  accurate  representation  of  relative  satellite  dynamics  in  a 
formation  flight  situation. 
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